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The existence of strongly bound excitons is one of the hallmarks of the newly discovered atomically 
thin semi-conductors. While it is understood that the large binding energy is mainly due to the 
weak dielectric screening in two dimensions (2D), a systematic investigation of the role of screening 
on 2D excitons is still lacking. Here we provide a critical assessment of a widely used 2D hydrogenic 
exciton model which assumes a dielectric function of the form e{q) = 1 + 27vaq, and we develop 
a quasi-2D model with a much broader applicability. Within the quasi-2D picture, electrons and 
holes are described as in-plane point charges with a finite extension in the perpendicular direction 
and their interaction is screened by a dielectric function with a non-linear ^-dependence which 
is computed ah-initio. The screened interaction is used in a generalized Mott-Wannier model to 
calculate exciton binding energies in both isolated and supported 2D materials. For isolated 2D 
materials, the quasi-2D treatment yields results almost identical to those of the strict 2D model and 
both are in good agreement with ah-initio many-body calculations. On the other hand, for more 
complex structures such as supported layers or layers embedded in a van der Waals heterostructure, 
the size of the exciton in reciprocal space extends well beyond the linear regime of the dielectric 
function and a quasi-2D description has to replace the 2D one. Our methodology has the merit of 
providing a seamless connection between the strict 2D limit of isolated monolayer materials and the 
more bulk-like screening characteristics of supported 2D materials or van der Waals heterostructures. 


I. INTRODUCTION 

Atomically thin semiconductors^ like graphene, hexag¬ 
onal boron-nitride (hBN), and M 0 S 2 are presently be¬ 
ing intensely studied due to their extraordinary opto¬ 
electronic properties. It is characteristic for these two- 
dimensional (2D) semiconductors that excitonic effects 
play a fundamental role, substantially modifying the op¬ 
tical spectrum by introducing states within the band gap 
that couple strongly to light and shift the onset of op¬ 
tical transitions to lower energies.Knowledge of the 
nature of the excitonic states is thus essential for device 
engineer ing^“^^. The well known Mott-Wannier model^^, 
which schematizes the exciton as a bound electron-hole 
pair interacting via a statically screened Coulomb interac¬ 
tion, is widely used to estimate exciton binding energies 
and radii in bulk materials. The main approximations 
behind the Mott-Wannier model are essentially three: 
(i) The real band structure is replaced by two parabolic 
bands, (ii) The microscopic shape of the conduction and 
valence band wave functions is neglected, (iii) The di¬ 
electric function is assumed to be local in real space, i.e 
g'-independent in reciprocal space. For 2D materials, the 
performance of the Mott-Wannier model and the valid¬ 
ity of the underlying approximations have still not been 
systematically investigated. The present work focuses on 
(iii), which is the only approximation where the role of 
the reduced dimensionality represents a qualitative dif¬ 
ference from the 3D case. 

For bulk semiconductors the macroscopic dielectric 
constant is defined as the limiting value of e(g) as g ^ 0. 
For a 2D semiconductor this definition cannot be straight¬ 


forwardly adopted since e(g = 0) = 1. In fact, for 2D sys¬ 
tems the dielectric function is strongly g-dependent and 
a more elaborate treatment of screening is required^^“^^. 
This issue has been treated by several authors^^’^^’^^’^^, 
who envisioned the 2D material as a strict 2D system, 
i.e. mathematically 2D, with a dielectric function of the 
form 

e 2 D(q) = 1 + 2TTaq, (1) 

where a is the 2D polarizability of the layer, which can 
be computed ab-initio. The screened electron-hole inter¬ 
action then follows 

O'TT 

W2D{cd = ( 2 ) 

where 27r/g is the 2D Fourier transform of 1/r. This form 
of interaction has the merit of leading to an analytical 
potential in real space and it has been successfully used 
to describe exciton binding energies and radii of several 
2D systems^^’^^’^^. We note that the form 1/g for the 
interaction and Eq. (1) for the dielectric function are 
consistent approximations which both become exact in 
the limit of vanishing thickness of the material, i.e. the 
strict 2D limit. However, to the best of our knowledge 
the validity range and limitations of these approximations 
have not previously been systematically explored. 

In this paper we relax the assumptions behind the 2D 
model adopting a microscopic approach that accounts for 
both the finite thickness of the layer and the full wave- 
vector dependence of the dielectric function. In the case 
of isolated monolayers our quasi-2D (Q2D) description 


2 




FIG. 1: Sketch of the (a) pure 2D and (b) quasi-2D 
Coulomb interaction. In the latter case the point 
charges can be thought as lines of charge extending 
along the thickness of the material. 


agrees well with the established strict 2D model provid¬ 
ing a justification for the latter. However, in the case 
of 2D layers supported on semi-infinite substrates or for 
thicker, i.e. few-layer, 2D materials, we find it important 
to account for the finite thickness and include the full 
non-linear g^-dependence of the dielectric function. In a 
recent paper we introduced a method for calculating the 
dielectric function of general layered materials (so-called 
van der Waals heterostructures^^“^^) where the dielectric 
functions of the individual layers are computed ab-initio 
and subsequently coupled together electrostatically^^. In 
the present work we use this method to compute the 
screened electron-hole interaction and solve the result¬ 
ing quasi-2D Mott-Wannier model for various types of 
heterostructures. We show that the exciton binding en¬ 
ergy and radius can be effectively tuned by controlling the 
screening via the heterostructure environment. Surpris¬ 
ingly we find that the transition from a strongly bound 
exciton in monolayer M 0 S 2 (binding energy of 0.6 eV) 
to a weakly bound exciton in bulk M 0 S 2 (binding energy 
of 0.15 eV) can be seamlessly described by the quasi-2D 
Mott-Wannier model accounting only for the change in 
the screening. 


II. THE QUASI-2D PICTURE 

Even though atomically thin semiconductors are re¬ 
ferred to as 2D materials they obviously do have a finite 
thickness. In this section, we show how the finite thick¬ 
ness can be accounted for within a 2D description. We 
shall refer to this description as the quasi-2D picture. 
To illustrate the concept, we consider the interaction be¬ 
tween two arbitrary charge distributions. 

(3) 

J |r-r'| 

In the case of two point charges confined to a 2D plane 


(see fig. 1(a)), each charge distribution is given by a delta 
function, i.e. pi{r\\) = qiS{T\\ — U,||), leading to an inter¬ 
action in reciprocal space: 

2'7r 

V2D(q||) = gi92|^- (4) 

Now we consider two charge distributions confined in 
a slab with finite thickness. We want to treat the real 
system, which is actually 3D, using an effective 2D de¬ 
scription. We do this by depicting the charge distribu¬ 
tions as lines of charge (fig. 1(b)). In other words, we 
assume that the charge densities are delta functions in¬ 
plane and have a certain distribution out-of-plane. The 
simplest approximation for the out-of-plane distribution 
is a step function of thickness d. This translates to 
p^(r||, z) = _ |^_ 2 ;q|)^ with zq the center of 

the material in the perpendicular direction, which leads 
to an interaction energy of the form (see appendix B): 


VQ2D(q||) 


47rgig2 

c^lqiiP 


1 - 

m\d 



(5) 

It is instructive to note that in the limit of q\\d 1 we 
recover the 2D potential, while for ^ 1 we get the 
Coulomb potential for a 3D system (calculated in-plane): 


VQ2D{q.\\) 


' 27rgig2 

i^iii 

47rgig2 

. |qiik 


q\\d < 1 
q\\d > 1 


( 6 ) 


III. SCREENED INTERACTION 


The (inverse) microscopic dielectric function gives the 
total potential to first order in the applied external po¬ 
tential, 

Vtot{r,uj) = J dr'e~^{r,r',Lo)Vext{r',uj), (7) 

here a harmonic time-dependence of the fields has been 
assumed. In standard ab-initio calculations for 3D pe¬ 
riodic systems, the dielectric matrix is calculated within 
the random phase approximation (RPA), which in plane 
waves representation takes the form: 

eGG'(q,w) = (5 gg' -^'(q + G)xGG'(q>‘^)> (8) 

with u(q + G) the Fourier transform of the Coulomb Po¬ 
tential and fhe non-interacting response function. For 
a 3D periodic system, the total potential resulting from 
a plane wave external potential has the form 

(9) 

where Fq(r, uj) is a lattice periodic function. Since usually 
we are interested in macroscopic fields, we define the 3D 
macroscopic dielectric function as 


1 

eM(q,w) 


yq(r,w) 


Vo 


eoo‘(q,w), 


( 10 ) 
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where {...)n denotes a spatial average over a unit cell. 
Note that in general eM(q, 7^ eoo(q, due to local 
field effects^^. 


A. Macroscopic Dielectric Function for 2D 
Semiconductors 


When eq. (10) is applied to an ab-initio calculation 
describing a 2D material as an infinite set of parallel 
sheets separated by a vacuum region of thickness L, 
eM(q, . This is a consequence of an av¬ 

eraging region much larger than the effective extension of 
the electron density around the material. The standard 
definition in eq. (10) becomes meaningless in this case, 
which is the reason why relatively different values for cm 
have been reported for monolayer M0S2 in the recent 
literature^’^^’^^. Therefore the definition of the macro¬ 
scopic dielectric function has to be revised accounting for 
the finite thickness. From the first equality in eq. (10), it 
is natural to substitute an average along the entire unit 
cell in the out of plane direction with an average over 
a confined region describing the actual extension of the 
electronic density. In practice, we average the in-plane 
coordinates (ry) over the unit cell area A and the out- 
of-plane coordinate from zq — d/2 to zq + d/2, where zq 
denotes the center of the material and d its width. The 
macroscopic dielectric function then becomes: 


Vq{r,Lo) 


A,d 


eQ2r>(q||,w) 


1^0 




G,zn sm{G±d/2) _-^ 


Gi. 


Gj_ o(*l|h‘^)> 

( 11 ) 


with eQQ,(q||,w) calculated from XQQ,(q||,w) according 
to the RPA expression in eq. (8). We stress that it is es¬ 
sential to use a truncated Coulomb potential in eq. (8) in 
order to decouple the layers in neighboring supercells^^. 
Note that we used the label Q2D since this definition 
of macroscopic dielectric function is consistent with the 
Q2D picture, as we show later on. As a rule of thumb we 
choose d to be the distance between the layers in the bulk 
form, but the results for excitons are not very sensitive 
to this choice as we shown in the next session. 

The q dependence of the static dielectric function is 
illustrated in fig. 2 for the case of monolayer hBN and 
M0S2. Without loss of generality, the q\\ values reported 
in the plot are chosen to be along the V — K direction. 
Indeed, further numerical tests show that the dielectric 
function is homogeneous, i.e. it is not significantly af¬ 
fected by different direction choices. In the low q\\ regime 
the dielectric function approaches one as expected for 2D 
materials^^. We mention in passing that the dielectric 
functions of a large collection of 2D materials is available 
in the Computational Materials Repository^^, see Refs. 
27 and 22 . 
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FIG. 2: Macroscopic dielectric functions for (a) hBN 
and (b) M0S2. The bulk(black), along with the Q2D 
(green) and 2D (blue) static dielectric functions are 
illustrated, the latter corresponding to the linear fits in 
the small q\\ region. For this calculations the q\\ values 
are taken along the V — K direction, but the 
homogeneity of the materials has been numerically 
verified. The parameters used in the linear response 
ab-initio calculation are discussed in section V. 


In the plots we also show the linear fit relevant for small 
q\\ as well as the bulk dielectric function. We see that for 
gyd ^ I a linear e is a viable approximation and we are 
in a 2D regime. In particular the 2D linear polarizability 
a can be calculated from the slope of the linear fit. On 
the other hand, when q\\d ^ I, the bulk behavior of the 
dielectric function is recovered. 


B. Screened Potential in Reciprocal Space 

To account for the screening in the charge-charge inter¬ 
action we modify eq. (3) introducing the dielectric func¬ 
tion: 


W 12 = [ drdr'dr' 

Jv 


„pi(r)e 7r,r")p2(r') 


( 12 ) 


In the following, we specialize to the case of electron-hole 
interaction. Assuming an in-plane delta function distri¬ 
bution and an unspecified ^(-dependence for the charge 
densities we can easily work out an expression for the 
screened interaction potential in reciprocal space: 


I poo 

^(‘l||) = |^ J dzdz'pe{z,ci\i)eoQ{z,z',ci\i)(l)h{z',ci\i). 

" (13) 

Here is the out-of-plane density distribution for 

the electron and 0/^(2), qy) is the out-of-plane potential 
generated by the hole. For details on how this poten¬ 
tial is calculated see appendix A. To study excitons in 
hBN and M0S2, we take the out-of-plane electron and 
hole distributions to be pe,h{z) = f^dr\\\'ipc,v K(i*y,^)P, 
with c and v the conduction and valence band indices 
respectively and K the high symmetry point of the first 
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FIG. 3: 2 ;-dependence of the total potentials (solid lines) coming from external perturbations (dashed lines) at 
different in-plane wave-vectors in the case of hBN ((a) and (c)) and M 0 S 2 ((b) and (d)). Left panels: the external 
perturbation is generated by a step function density distribution (insets). Right panels: the external perturbation is 
generated by the actual hole out-of-plane density distribution (insets), which is calculated as 
Ph{z) = X 4 ^^11 l'^vK(r||, z)\‘^ with V indicating the valence band and K the high symmetry point of the first Brillouin 

zone. In all cases the density distributions are normalized to 1. 


Brillouin zone, since for both materials that is where the 
lowest bound exciton is localized^Furthermore, in 
eq. (13) we have introduced a mixed representation for 
the dielectric function, specifically: 


G±G'^ 

(14) 

Note that taking G\\ = Gj| =0 corresponds to an in-plane 
macroscopic dielectric function, which also accounts for 
local field effects. 


general expression eq. (13) reduces to (see appendix C): 


Wq2d{^\\) 


-1 / N 

rf|q|||2eQ2D(q||)x 

L m\d 

= eQ2D(q||)^Q2r>(q||), 



(15) 


where ^Q 2 Di^\\) ^be static version of the macroscopic 
dielectric function defined in eq. (11). We thus see that 
^Q 2 D is the natural dielectric function to be used in the 
quasi-2D picture. 


To illustrate the effect of screening, fig. 3 shows how 
a potential generated by either the step function den¬ 
sity distribution or the actual hole density distribution 
is screened by hBN and MoS2. In all cases the density 
distribution is normalized to 1. The possibility of using 
either the actual electron/hole out-of-plane density dis¬ 
tribution (fig. 4) or a simply step-function gives us two 
different approximations to calculate the screened inter¬ 
action within the Q2D picture. 

In the case of step-function density distributions, we 
can find an analytic expression for the screened potential 
in eq. (13), if we make a further approximation. Indeed, if 
instead of considering the full ^-dependence of qy) 
we take its average value within a region of thickness d 
around the layer, and then screen the resulting constant 
potential by the full ^(-dependent dielectric function, the 




^(A) ^(A) 

FIG. 4: Valence (red) and conduction (green) band 
densities for (a) hBN and (b) M 0 S 2 calculated at the K 

point. 
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FIG. 5 : Macroscopic dielectric functions for (a) hBN 
and (b) M0S2. The different dielectric functions are 
calculated with the three different approaches explained 
in the text: dielectric function from actual electron and 
hole distributions (magenta), dielectric function from 
step function distributions (cyan) and Q2D dielectric 
function (green). The vertical line represent the radius 
of the exciton in reciprocal space. 

For each of the two different Q2D models for the 
screened electron-hole interaction, we can associate a 
Q2D dielectric function defined as the ratio between the 
bare and the screened potential: 

,7 w (p?(qii)l</>Z(qii)) 

" (pI(qii)|eoo^(^,5',qii)l'/>yqii))’ 

where for simplicity we have used a bracket notation for 
the integration over ^ and 7 = steps,wfs indicates whether 
the potentials are calculated from step functions or actual 
electron and hole density distributions. Figure 5 shows a 
comparison of the two dielectric functions thus obtained 
together with eg2D from eq. (11) for hBN and M0S2. 
Clearly the curves perfectly agree in the low q\\ regime, 
while deviations appear for higher values. This observa¬ 
tion is consistent with the fact that for small wave vectors 
the total potentials are flat, and therefore well approxi¬ 
mated by the Q2D average value (see fig. 3 ). As we show 
later, the relevant q\\ region for the screening is the one be¬ 
low the the black vertical line representing the inverse ex¬ 
citon radius, calculated from the ab-initio Bethe-Salpeter 
Equation (BSE) (see section V). Therefore the three dif¬ 
ferent Q2D approaches can be considered equivalent when 
dealing with excitons in these monolayer materials. 

C. Screened potential in real space 

To obtain the form of the screened potentials in real 
space we Eourier transform eq. ( 13 ): 

IT. / ^ 2 /■“ Mq\rii\) _1 

Wq2d{v\\) =-^ dq --- eQ2Diq)x 

dJo q 

i_ 2 w, 

[ qd V 2 yj ’ 


EIG. 6: Screened Q2D and 2D potentials for (a) hBN 
and (b) M0S2. The potentials are calculated 
numerically starting from the macroscopic dielectric 
functions in fig. 2 and using eq. ( 17 ) and eq. ( 18 ) 
respectively. The bare Q2D curves are calculated using 
the first equation but neglecting the screening. 


where Jo{x) is the zeroth order Bessel function and where 
we used that the dielectric function is isotropic. This is 
the quasi- 2 D potential which can be compared to its strict 
2 D counterpart defined in eq. ( 2 )^^: 

W 2 D{rii) = ^ [Hoix) - , ( 18 ) 

where Ho{x) and No{x) are the Struve and Neumann 
functions respectively. We stress here that the parameter 
a can be estimated from the slope of the fit in fig. 2. We 
note that while this procedure of calculating the 2 D po¬ 
larizability differs from the standard one, it is equivalent. 
In the case of M0S2, for example, we obtain a value of 
5 . 9 A which agrees well with the value of 6.6A obtained 
in the literature^^. 

In fig. 6 we report the numerical results for different 
potentials: the bare Q2D (black) obtained from eq. ( 17 ) 
setting eQ2D to 1 , the screened Q2D (green) obtained 
from the same equation but including the screening as 
^Q2D and the screened 2 D (blue) calculated from eq. ( 18 ). 
The results are shown for both hBN and M0S2. 

We note that the bare Q2D interaction agrees with 
— 1/r beyond a distance given by the layer thickness d. 
Eurthermore we see that increasing the layer thickness 
(going from hBN to M0S2) reduces the bare Q2D inter¬ 
action strength as expected from eq. ( 17 ). Including the 
screening reduces the interaction strength even further. 
The reduction is most significant when using the linear 
dielectric function (strict 2 D screening) as expected from 
fig. 2 which shows that 621) (^) > ^Q2d{q) for all q. We 
see that, apart from a significant deviation for electron- 
hole separation smaller than roughly lA , the 2D and 
Q2D screened potentials agree and both show a logarith¬ 
mic dependence for r ^ 0 . Eor distances larger than the 
layer thickness, all the potentials (screened and bare) ap¬ 
proach the same value (—1/r), meaning that screening is 
completely absent in the asymptotic limit. 
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FIG. 7: Variation of the macroscopic dielectric function 
and effective potential in M 0 S 2 due to the change in the 
thickness d of the averaging region in the Q2D model. 

The continuous black lines are relative to 
d = 6.29A (the interlayer distance in the bulk), while 
the dashed lines delimiting the shaded region are 
calculated with a variation of ±10% in d. 


D. Importance of the Thickness Parameter 

We now return to the problem of choosing the exter¬ 
nal parameter d entering the Q2D dielectric function. In 
fig. 7 we show the Q2D dielectric function and the corre¬ 
sponding potentials when d is varied by ±10% around the 
interlayer distance in bulk M 0 S 2 . To the left of the maxi¬ 
mum, eQ 2 D is insensitive to d since the induced potential 
is constant over the averaging region. Also in the high 
limit, eQ 2 D is not affected. This is because for these wave 
vectors the induced potential is in practice negligible. In 
general, increasing (decreasing) d, decreases (increases) 
eQ 2 D in the large wave vector region. Despite the fact 
that the variation in the dielectric function is fairly vis¬ 
ible for intermediate g-values, the screened potential is 
barely modified. This is because the bare Q2D potential 
shows an opposite dependence on d, such that the product 
WQ 2 D{q) = eQ 2 D( 9 )^Q 2 D( 9 ) stays essentially unchanged. 
In terms of exciton binding energies we have found that 
a ±10% variation in d leads to a correction of less than 
O.OleV. 


IV. GENERALIZED MOTT-WANNIER MODEL 

An accurate description of excitonic effects requires 
the solution of a computationally demanding many-body 
problem, namely the Bethe-Salpeter equation (BSE)^^’^^. 
However, it is well known for 3D systems that a satisfy¬ 
ing qualitative description can be obtained modelling the 
exciton as a hydrogenic atom constituted by an excited 
electron-hole pair interacting via a statically screened 
Coulomb interaction. In this section we generalize such 
a model to the Q2D case. 

The Bethe-Salpeter two particle Hamiltonian for a 2D 


periodic system is given by: 


(^n2ki+q|| ^niki )^nin3^n2n4^kik2 ± 
n3n4k2 (19) 

~ (/n2ki+q|| ~ /niki )A"nin2ki (qy )? 

n3n4k2 


where rii are band indices, vectors in the first 2D Bril- 
louin zone and qy is the in-plane momentum transfer, or 
exciton center-of-mass momentum. In the following we 
specialise to the case of vertical transitions, i.e. qy = 0. 
K is the kernel containing the exchange and the screened 
direct Coulomb interaction. This Hamiltonian describes 
scattering processes between two electron-hole pairs ex¬ 
cited by an external perturbation. In general these pro¬ 
cesses should involve all the occupied and unoccupied 
states in the spectrum; however, when the conduction 
and valence bands are well separated from the remain¬ 
ing bands, it is often a good approximation to include 
only the valence and conduction band states. Together 
with the Tanm-Dancoff approximation this assumption 
allows us to express the resonant part of the two-particle 
Hamiltonian as: 


= (eck' - e.k)4k' - i^«ck . 
vck' ^ck' 


The kernel is given by 


( 20 ) 


vck' J A 

2 J drdr''tpy]^{r)i;*i^{r)v{r,r')'tp*^,{r')'tpck'{r') 

( 21 ) 

where IV^ak)? with a = {v,c)^ represents single parti¬ 
cle Bloch states for the valence and conduction band, 
IT(r,r') = f is the screened interaction and 

i;(r, r') = is the bare Coulomb interaction. The 

first term on the right hand side of eq. (21) is the direct 
screened electron-hole interaction while the second is the 
Coulomb exchange. Our full ab-initio solution of the BSE 
shows that the exchange term only slightly decreases the 
exciton binding energy by 0.08 eV and 0.02 eV for hBN 
and M 0 S 2 , respectively. This amounts to less than 5% 
of the total binding energy and we therefore neglect the 
exchange contribution in the rest of the paper. 

Throughout the BZ we consider the valence and con¬ 
duction band wave functions to be plane waves in the 

in-plane direction and in the out-of-plane direction equal 

1 /2 

to dry |'0Q,x(i*y, 2 ))p) up to a normaliza¬ 

tion factor and with a = v^c. With this approximation 
and proceeding as for eq. (13), the interaction matrix be¬ 
comes: 


if.ck =4^(|k-k'|), (22) 

where IF(|k|) is the screened potential in eq. (13) which 
can be evaluated in the various ways described in the pre¬ 
vious section, depending on the level of approximation. 
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TABLE I: Geometry and effective masses. 


Material 

a (A) L(A) d (A) 

(a.u.) 

M 0 S 2 

3.20 23.0 

6.29 

0.27 

hBN 

2.50 23.0 

3.22 

0.37 


Completely analogue to the 3D case we can intro¬ 
duce the envelope function E(r||), defined as E(r||) = 
A(k), with A(k) excitonic weights in recipro¬ 
cal space, and arrive at an eigenvalue problem of a 2D 
hydrogenic atom, 



r 

^ 2D 


+ W{^) 


F(r||) = EbF{rii), 


(23) 


where flex is the exciton effective mass, calculated from 
the hole and electron masses according to = m~^ -h 


FIG. 8: Convergence plot for the binding energy 
obtained from the BSE solution against the k-points 
density. Extrapolation to infinite k-point sampling is 
shown. In the BSE the exchange contribution is left out 
The horizontal dashed lines show the results given by 
the Q2D model. 


V. EXCITON BINDING ENERGIES OF 
ISOLATED MONOLAYERS 

In this section we investigate the performance of the 
Mott-Wannier model in eq. (23) for the calculation of 
binding energies of the lowest bound exciton in hBN and 
M 0 S 2 . 

A. Ab-initio Calculation Details 

In order to solve eq. (23) with either the Q2D or 2D 
potentials, we need to calculate the dielectric matrix. 
We describe the two materials with a supercell technique 
and we optimize the structure using the LDA exchange- 
correlation potential; geometrical details are provided in 
table 1. To calculate the non-interacting response func¬ 
tion we use 150eE cut-off energy for the reciprocal lattice 
vectors G and G' in order to account for local field ef¬ 
fects. We construct from local density approximation 
(LDA) wave functions and energies, and we then get the 
dielectric matrix using a truncated Goulomb potential in 
order to avoid interaction between supercells^^. The di¬ 
electric matrix is calculated on a 60 x 60 k-points grid. 
Since it turns out that the exciton binding energy is sen¬ 
sitive to the low wave-vector behavior of the screening, 
we use an expansion of the density response function 
around q\\ = 0 in order to calculate the dielectric ma¬ 
trix in the small limit. All calculations are performed 
with the GPAW code^^’^^, which is based on the projector 
augmented wave method. Details on the implementation 
of the linear response code can be found in Ref. 34. We 
mention that the dielectric functions of more than fifty 
2D materials calculated in this fashion are available in 
the GMR^^. The exciton masses as computed from the 
LDA band structure are given in table 1. 

To obtain the lowest bound exciton we numerically 
solve the Mott-Wannier equation on a logarithmic grid. 


With this method we are able to converge the lowest 
eigenvalue with a precision of 0.002eV. For benchmark, 
we perform BSE calculations using the GPAW code. For 
the screening of the electron-hole interaction we use the 
static dielectric function evaluated with the same param¬ 
eters employed in the linear response calculation. The 
particle-hole states of the BSE Hamiltonian are con¬ 
structed from a single LDA valence and conduction band. 
To compare directly to our model, all the BSE calcula¬ 
tions are performed neglecting the exchange part of the 
kernel. We stress that the binding energy of the first exci¬ 
ton changes by less than O.OleE if the BSE Hamiltonian 
is constructed from the four highest and four lowest con¬ 
duction bands. As reported previously^BSE binding 
energies in 2D materials are extremely sensitive to the k- 
point grid. We therefore perform BSE calculations with 
up to 60 X 60 /c-points, for which we get binding energies 
of 2.07eV and 0.54eV for hBN and M 0 S 2 respectively. 
Furthermore, assuming a linear dependence of the bind¬ 
ing energy with respect to 1/pkpts^ we extrapolate the 
results to infinite /c-points sampling (see fig. 8). 


B. Results 

The values for binding energy of the lowest bound ex¬ 
citon obtained with the different models for the screened 


TABLE H: Numerical values for energy of the lowest 
bound excitonic state at the direct gap. Both the BSE 
and the models are based on LDA ab-initio calculations. 
The exchange contribution is not included. 

Eg^^(eV) Eg^^(eV) Eg^feV) E^^^(eV) 

hBN 2.05 2.35 2.34 2.23 2.29 

MoS2 0.43 0.61 0.60 0.57 0.59 



















electron-hole interaction along with the extrapolation 
from the BSE are reported in table 11. We first observe 
that there is practically no difference in the binding ener¬ 
gies obtained from the Mott-Wannier model using either 
the Q2D or 2D screened potentials. Moreover, the re¬ 
sult from the Mott-Wannier model(s) are within 0.3eV 
and O.lSeV of the BSE result for hBN and M 0 S 2 respec¬ 
tively. We consider this a reasonable agreement given the 
simplicity of the model. 

In table II we also report the binding energies obtained 
when the electron-hole interaction is calculated numeri¬ 
cally from eq. (13) using step-functions and actual elec¬ 
tron and hole density distributions. As pointed out in the 
discussion of fig. 5 we expect these two other approaches 
to give the same description of excitons. Indeed, the bind¬ 
ing energies we obtained are in perfect agreement with the 
Q2D and 2D model results. 

The agreement between the Q2D and 2D descriptions 
can be understood by looking at the g-space extension of 
the lowest bound exciton wave function shown in fig. 9. 
We see that for both hBN and M 0 S 2 the exciton is con¬ 
fined to a rather narrow region around the K-point. A 
localization of the exciton in g-space means that the rel¬ 
evant contribution to the electron-hole interaction comes 
from the low wave-vector regime. Erom the calculated 
excitonic wavefunctions in real space we obtain inverse 
exciton radii of 0.29A“^ for hBN and O.OTA”^ for M 0 S 2 . 
Both of these values are comparable to 1/d (0.3lA“^ and 
0.16A“^, respectively). As we have seen previously, in 
this limit the Q2D screened potential reduces to the strict 
2D potential explaining the similarity of the binding en¬ 
ergies obtained with the two descriptions. 

To conclude this section, we notice that in the eval¬ 
uation of the screened electron-hole interaction, we ne¬ 
glected the in-plane spatial variation of the conduction 
and valence band wavefunctions. The validity of this ap¬ 
proximation can be checked by performing a BSE calcu¬ 
lation where the screened interaction is evaluated using 
a dielectric matrix Cqq, where all matrix elements ex¬ 
cept for those where Gy = Gj| =0 are set to zero. In 
other words, we neglect all the in-plane high frequency 
spatial variations of the wave functions. With this con¬ 
striction we obtain a binding energy of 2.21eE for hBN 
and 0.44 for M 0 S 2 . The neglect of in-plane variations of 
the wave functions is thus responsible for 0.15 eV (hBN) 
and 0.01 eV (M 0 S 2 ) of the observed discrepancy between 
the Mott-Wannier model and the full BSE calculation. 


VI. EXCITONS IN LAYERED STRUCTURES 

In this section, we show that a linear approximation for 
the dielectric function breaks down when applied to ex¬ 
citons in multi-layered structures. While it is possible to 
include the non-linear q-dependence of the dielectric func¬ 
tion within a strict 2D model, the Q2D description turns 
out to be necessary to quantitatively capture screening 
effects. 


A. The Quantum Electrostatic Heterostructure 
(QEH) Model 

In order to calculate exciton binding energies in a lay¬ 
ered structure we first need the dielectric function. This 
can be obtained using the quantum-electrostatic het¬ 
erostructure (QEH) model that we introduced recently^^. 
In brief, the underlying procedure in the calculation of the 
dielectric function can be divided in two parts. In the first 
part the full density response function of each isolated 
layer, calculated from first principles, is used to obtain 
the monopole/dipole components of the density response 
function as well as the spatial profile of the electron den¬ 
sities in the ^-direction induced by a monopole/dipole 
field. We refer to these data sets as the dielectric build¬ 
ing block of the individual layer. In the second part, the 
dielectric building blocks are coupled together via the 
Coulomb interaction to give the dielectric function for 
the full structure. The dielectric matrix obtained from 
the QEH model can be used to obtain the electron-hole 
interaction according to 

W{q\\) = Pl{q\\) r^(9||) (24) 

where p ^ electron density (hole induced poten¬ 

tial) vector expressed in a basis set of monopole/dipole 
densities (potentials). The basis set of induced densi¬ 
ties and potentials is also used as (left and right) ba¬ 
sis functions for representing To be more explicit 

an arbitrary density vector p can be represented as 
^ = [piM, piD, p 2 M, p 2 D, • • • , PnM, Pud] where pia, with 
a = is the induced monopole/dipole density at 

the layer i. A completely equivalent expression can be 
formulated for the induced potentials. 

It is clear that the equation above is just a sim¬ 
ple rewriting of eq. (13) in terms of a minimal 
monopole/dipole basis. We point out that this formalism 
takes the finite extension of the layers in the out-of-plane 
direction into account, and is therefore consistent with 
the Q2D picture described in the previous sections. In 
Ref. 22 we showed, based on the comparison with full 
ab-initio calculations, that the monopole/dipole basis is 
sufficient to obtaining an accurate description of the di¬ 
electric and plasmonic properties of different layered het¬ 
erostructures. 


B. Breakdown of the Linear Screening Model 

As an example we consider two different types of het¬ 
erostructures. The first, which we refer to as “on-top”, 
consists of M 0 S 2 on top of n layers of hBN. The second, 
which we refer to as “sandwich”, consists of an M 0 S 2 
layer encapsulated in n layers of hBN, see fig. 10 panel 
(a) and (c). The interlayer distance between M 0 S 2 and 
hBN is set to 5.1 A. In fig. 10 panel (b) and (d) we show 
the dielectric function of the M 0 S 2 layer as well as the 
linear approximation as a function of the in-plane wave 
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FIG. 9: LDA band structure and exciton weights for (a) hBN and (b) M 0 S 2 . In both materials the exciton is well 
localized at the K point. The excitonic weights are calculated as the absolute value squared of the eigenvector of the 
two particle BSE Hamiltonian associated to the lowest bound exciton. In red the parabolic bands used in the 
Mott-Wannier model. The values for the electron and hole masses are 0.93a.i4. and 0.62a.i4. for hBN and O.Gla.rt. 

and 0.49a.14. for M 0 S 2 . 


vector for different number of hBN layers. The effective 
dielectric function of M 0 S 2 in the heterostructure is de- 
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FIG. 10: Left panels: The on-top (a) and sandwich (c) 
arrangements of the MoS 2 /hBN heterostructures. Right 
panels: effective dielectric function (full line) for the 
on-top (b) and sandwich (d) configurations. The linear 
approximations to the dielectric function is shown by 
dashed lines. The shaded regions in (b) and (d) 
represent the range of inverse exciton radii found for the 
considered structures. The values below these regions 
are relevant for screening the electron-hole interaction 
and for the thicker structures this region extends 
beyond the linear regime of e{q). 


fined along the lines of eq. (16): 


e(9||) = 


£l(g||) ^fe(g||) 

p\{<i\\) g~\<i\\) 


(25) 


which gives the ratio of the bare to the screened interac¬ 
tion between an electron and a hole in the M 0 S 2 layer. 

From fig. 10, we notice that adding hBN layers to M 0 S 2 
changes the shape of the dielectric function introducing a 
pronounced feature that shifts towards low q\\ as the num¬ 
ber of hBN layers is increased. This shoulder-like feature 
can be explained as an interplay between the 3D and 2D 
screening character. When more hBN layers are added to 
the heterostructure, the system tends toward a bulk limit, 
where the dielectric function is larger than 1 for q\\ = 0. 
However, the heterostructure has a finite thickness d and, 
as required by the 2D limit, when gy ^ 1/d the dielec¬ 
tric function is 2D like and becomes 1 for gy = 0. This 
leads to a sharp drop in the dielectric function, which be¬ 
comes steeper as the thickness of the heterostructure is 
increased, explaining the appearance of the shoulder-like 
feature. It is clear, from fig. 10, that the main change 
in the dielectric function is caused by the nearest layers 
of hBN. Adding more layers only causes a slight varia¬ 
tion. Obviously, this is due to the fact that hBN is less 
effective at screening the electron-hole interaction as the 
distance from M 0 S 2 is increased. For the same reason, 
the screening is more pronounced in the sandwich config¬ 
uration than in on-top configuration. 

We then proceed to calculate the binding energy of the 
lowest bound exciton in the M 0 S 2 layer for the two differ¬ 
ent configurations using both the full wave vector depen¬ 
dent dielectric function (quasi-2D) and its linear approxi¬ 
mation (strict 2D). The results are shown in fig. 11 panel 
(a) and (b). When the full dielectric function is used. 
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FIG. 11: Energy and radius of the lowest bound exciton 
for the (a,c) on-top and (b,d) sandwich configuration as 
function of the number of hBN layers obtained from the 
Q2D (green) and 2D (blue) approaches. 


the binding energy converges towards 0.40eE and 0.31eE 
for the on-top and sandwich configurations, respectively. 
These values represent the bulk limits, i.e. M 0 S 2 on a 
hBN substrate and M 0 S 2 encapsulated by two seminfi- 
nite stacks of hBN. The reduction in binding energy of 
0.2 eV for the on-top configuration is in good agreement 
with the experimentally determined change in exciton en¬ 
ergy of WS 2 when adsorbed on Si 02 ^^ (the bulk dielec¬ 
tric constants of Si02 and hBN are similar). In contrast, 
the assumption of linear dielectric screening completely 
fails in estimating the exciton binding energies. Indeed, 
it quickly diverges from the Q2D results yielding much 
too small binding energies. This behaviour results from 
the continuously increasing slope of the dielectric function 
eventually arriving at a condition of perfect screening (in¬ 
finite slope). 

Figure 11 panels (c) and (d) show the exciton radii 
obtained from the Q2D and 2D models. Interestingly, 
for the Q2D model the increase in the exciton radius 
due to the screening from the hBN is only 10% and 30% 
for the on-top and sandwich configurations, respectively. 
The range of the inverse exciton radius is indicated by a 
shaded region in fig. 10 panel(c) and (d). As we demon¬ 
strated in the previous section, the relevant gy for the 
screening lie mainly below the inverse exciton radius. In¬ 
spection of fig. 10 clearly indicates that in this regime 
the linear approximation overshoots the full wave vec¬ 
tor dependent dielectric function and it gets worse as the 
number of layers is increased. 


C. Limitations of the 2D picture in Layered 
Structures 

In the previous paragraph we showed that the assump¬ 
tion of linear screening, i.e. eq. (1), breaks down when 
the screening from the environment is included. It is, of 
course, possible within the 2D picture to couple a stack of 
2D materials, each described by a linear dielectric func¬ 
tion, using the QEH model. In this section we explore 
the validity of such an approach using the Q2D results 
obtained in the previous section as a reference. 

We model the layered structure as infinitesimally thin 
planes described by 2D building blocks, as opposed to the 
Q2D ones used previously, and couple them electrostati¬ 
cally via the QEH. While it is straightforward to define 
multipole response function and induced density compo¬ 
nents when a finite thickness is considered, in 2D only the 
monopole components have an obvious definition. Within 
the 2D picture, the monopole induced density is described 
by a delta function centered at the layer position Zi. The 
component of the 2D response function of the (isolated) 
layer may be obtained from the corresponding 2D dielec¬ 
tric function in eq. (1): 

xlhiQw) = Ur [e2“2(9||) - 1 ] 


1 + 27raq\\ 

For strict 2D layers, the Coulomb interaction between 
monopole charge densities in layers at Zi and z^ takes the 
form 

ViM,kM{q\\) = --, (27) 

which reduces to the standard 2D Coulomb potential in 
reciprocal space for coupling within the layer. 

To test the QEH with 2D building blocks, we consider 
the “on-top” structure of the previous paragraph and in 
fig. 12 we report the effective dielectric function and en¬ 
ergy of the lowest bound exciton as a function of the num¬ 
ber of hBN layers. It is clear that the 2D dielectric func¬ 
tion of the supported M 0 S 2 deviates significantly from 
the Q2D result (see fig. 10 panel (b)) for larger gy. How¬ 
ever, for smaller gy the 2D function actually reproduces 
qualitatively the non-linear structure of the Q2D result. 
In terms of exciton binding energy, we observe a conver¬ 
gence to a finite value when the number of layers of hBN is 
increased, but the reduction in binding energy compared 
to the free-standing layer is 50% smaller than the reduc¬ 
tion obtained with the Q2D approach. The underestima¬ 
tion of the screening can be ascribed essentially to two 
reasons. First, the potential generated by a 2D induced 
density decays faster than the actual one, making the 
neighboring layers less effective at screening the electron- 
hole interaction. Second, the dipole response of the lay¬ 
ers, which would increase the screening even more, is not 
included. In particular we mention that within the Q2D 
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FIG. 12: (a) Effective dielectric function (full line) and 
(b) energy of the lowest bound exciton for the on-top 
MoS 2 -hBN configuration, calculated with the QEH 
model based on a 2D description of the layers. The 
linear approximations to the dielectric function is shown 
by dashed lines in panel (a), along with the range of 
inverse exciton radii found for the considered structures 
represented by the shaded region. 
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EIG. 13: Energy of the lowest bound exciton for M 0 S 2 
incapsulated in M 0 S 2 layers as function of the total 
number of M 0 S 2 layers obtained from the Q2D (green) 
and 2D (blue) approaches. With the Q2D approach we 
can clearly see the transition from the mono-layer 
exciton to bulk one. 


approach, removing the dipole contribution increases the 
converged value of the binding energy by O.OTeV. To con¬ 
clude this paragraph, we have shown that even though 
the 2D model does capture the essential non-linear shape 
of e(g||) in the small q\\ regime, it underestimates the ef¬ 
fect of environmental screening and consequently predicts 
too small changes in exciton binding energies due to sub¬ 
strate effects. 


D. Transition from 2D to 3D Excitons in M 0 S 2 


As a final example, we study the 2D to 3D transition 
of the exciton in M 0 S 2 . In layered bulk materials, the 
Mott-Wannier equation can be written as follow: 




+ lF(r) 


F(r) = E,F{v), 


(28) 


where typically the exciton mass in the out of plane 
direction is much higher than the in-plane directions 
^ Consequently, we can neglect the out-of- 

plane component of the kinetic energy and be left with the 
2D Mott-Wannier model. Additionally, the in-plane ef¬ 
fective mass does not vary considerably going from mono- 
layer to bulk M 0 S 2 as shown in Ref. 36. Therefore, 
the main difference between the physics of excitons in 
monolayer and layered bulk is contained in the screened 
potential rather than the geometric confinement. Based 
on this, it is tempting to model the bulk exciton as an 
electron-hole pair confined to a single layer but with an 
interaction screened by the bulk environment. To test 
this we consider a multi-layer M 0 S 2 structure and calcu¬ 
late the binding energy of an exciton localized in the cen¬ 
tral M 0 S 2 layer using the Q2D Mott-Wannier model with 
screened potential calculated from the QEH model. The 
results for the binding energy as a function of the number 
of M 0 S 2 layers are plotted in fig. 13. As expected, the 


reduction of the exciton binding energy is larger when 
the monolayer is embedded in M 0 S 2 than in the case of 
hBN (fig. 11 panel(b)). Amazingly, the binding energy 
converges towards a value of 0.16eE only 0.03eE higher 
than previously reported ab-initio value for bull MoS 2 ^^. 
This shows that the different nature of excitons in 2D 
and layered 3D materials is mainly caused by the screen¬ 
ing while quantum confinement plays a minor role. 


VII. CONCLUSIONS 


In this work we have presented a systematic study of 
the screening properties of two-dimensional semiconduc¬ 
tors and layered structures. Taking into account the finite 
extension of the 2D material in the out-of-plane direction, 
we have proposed a general quasi-2D picture to describe 
the screened electron-hole interaction in the context of 
excitons. We have shown that, in the case of isolated 
layers, the excitons are typically large enough that the 
screening can be described by a linear dielectric function 
consistent with a strict 2D picture. On the other hand, 
for multi-layer structures where the screening properties 
are intermediate between the 2D and 3D regimes, it is 
essential to include the non-linear g^-dependence of the 
dielectric function. If this is done and a quasi-2D descrip¬ 
tion is employed, very satisfactory results are obtained for 
both monolayer and multi-layer structures using the same 
theoretical framework. In combination with a recently 
introduced scheme for computing dielectric functions of 
layered materials^^, this makes it possible to model exci¬ 
ton physics in general van der Waals heterostructures at 
very low computational cost. 
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Appendix A: Poisson’s equation for lines of charge 

Charges in 2D materials can be depicted as lines ex¬ 
tending over the thickness of the layer. The potential 
generated by a line of charge can be calculated from the 
Poisson equation: 


V^(^(r) = —47rp(r). (Al) 


Because of the cylindrical symmetry of the line of charge, 
it is convenient to Fourier transform in the in-plane di¬ 
rection and rewrite eq. (Al) as: 



^’(qih^:) = -47rp(q||,2:). 


(A2) 


For a line of charge, the density distribution can be 
separated as an in-plane delta function and an out-of- 
plane function p{z) and therefore its in-plane Fourier 
transform would read p(q||, 2 :) = p{z). From the 

structure of eq. (A2) and the form of the Fourier trans¬ 
formed density, it is convenient to write the potential 

as = q^jy!-(/)(z,q||). Note that is the 

Fourier transformed solution for the Poisson equation for 
a point charge in a 2D plane, therefore we can consider 
(j){z^Q[_\\) as the out-of-plane component of the potential. 
Plugging ip{z^ci\\) and p(q||, 2 :) in eq. (A2), we finally ob¬ 
tain the Poisson equation for the out-of-plane potential 
generated by a line of charge: 

^2 

- |q|y^!.(2:,q||) = -47r|q|||(A3) 


Appendix B: Unscreened Q2D Interaction 

In this appendix we derive the expression for the Q2D 
unscreened charge-charge interaction in eq. (5). Accord¬ 
ing to our Q2D picture and assuming a charge distribu¬ 
tion pi(Y\\,z) = — jz —zol), the unscreened 

charge-charge interaction in reciprocal space can be writ¬ 
ten as: 

T. . ^ f , e{i-\z-z^\)e^rn 

VQ2D{^.\\) =qiq2 - 2 - 

|r-r'| d 

(Bl) 



To proceed we notice that the integral in the second line 
can be interpreted as the potential generated by an in¬ 
plane Fourier transformed charge distribution p(q||, z) = 

£( 2 — 1^ ^ol) g-^q||-r|| analytic form can be obtained 

solving eq. (A2) as illustrated in appendix A: 


, , 47re-^^r^ii 

fB2') 

fl - cosh(|q|| Ijz' - Zo\) \z' - Zo\ < ^ 

|e-|qiill^-^olsinh(|q|||d/2) > f ‘ 

Plugging this result in eq. (Bl) and integrating in-plane 
and along z separately, we recover the expression eq. (5) 


Appendix C: Screened Q2D Interaction 


In the following we show how to derive the expression 
for the Q2D screened electron-hole interaction reported in 
eq. (15). For charge distributions of the kind pi{T\\,z) = 

— jz — zol), the screened interaction reads: 


WQ2D(q||) = - / drdr 


i 


^,d{l-\z-zo\)e^^r- 


Jv 


dr" 


9(1 - \z" 


zo\)e~ 


mil -r 


(Cl) 


As done in appendix B we can interpret the integral in the 
second line as the potential in eq. (B2). In order to keep 
the calculation analytic we approximate (p{q _\\, z) with its 
average inside the slab in the out-of-plane direction, in 
formula: 




ZQ-\-dj2 


ZQ—dl2 

47re-'qii-^ii 


dzip{q_\\,z) 


d2|q|||2 

g-myrii 


1 - 


h\\d 


d 


VQ2_D(q||) 


e-|qiiM/2 sinh(q||d/2) 


(C2) 


Inserting the last expression in eq. (C2) and integrating 
in-plane we get: 


2 nZo+dl2 nZo-\-l^l'Z 

Wq2d{(1\\) = Vq2d{cii\)- dz f^ooi^^z') 

^ Jzo-dl2 Jzn-Ll2 


zo-\-L/2 


ZQ-LI2 


= l^Q2D(q||)eQ2D(q||)- 

(C3) 

















13 


^ Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Cole¬ 
man, and M. S. Strano, “Electronics and optoelectronics 
of two-dimensional transition metal dichalcogenides,” Na¬ 
ture nanotechnology^ voL 7, no. 11, pp. 699-712, 2012. 

^ K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, 
“Atomically thin M 0 S 2 : A new direct-gap semiconductor,” 
Phys. Rev. Lett., vol. 105, p. 136805, Sep 2010. 

^ A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. 
Chim, G. Galli, and F. Wang, “Emerging photolumines¬ 
cence in monolayer M 0 S 2 ,” Nano letters, vol. 10, no. 4, 
pp. 1271-1275, 2010. 

^ A. Ramasubramaniam, “Large excitonic effects in mono- 
layers of molybdenum and tungsten dichalcogenides,” 
Phys. Rev. B, vol. 86, p. 115409, Sep 2012. 

^ D. Y. Qiu, F. H. da Jornada, and S. G. Louie, “Optical 
spectrum of M 0 S 2 : Many-body effects and diversity of ex- 
citon states,” Phys. Rev. Lett., vol. Ill, p. 216805, Nov 

2013. 

® M. M. Ugeda, A. J. Bradley, S.-F. Shi, H. Felipe, Y. Zhang, 
D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, 
et al., “Giant bandgap renormalization and excitonic ef¬ 
fects in a monolayer transition metal dichalcogenide semi¬ 
conductor,” Nature materials, 2014. 

^ K. He, N. Kumar, L. Zhao, Z. Wang, K. F. Mak, H. Zhao, 
and J. Shan, “Tightly bound excitons in monolayer WSe 2 ” 
Phys. Rev. Lett., vol. 113, p. 026803, Jul 2014. 

® D. Jariwala, V. K. Sangwan, L. J. Lauhon, T. J. Marks, and 
M. C. Hersam, “Emerging device applications for semicon¬ 
ducting two-dimensional transition metal dichalcogenides,” 
ACS nano, vol. 8, no. 2, pp. 1102-1120, 2014. 

^ M. Bernard!, M. Palummo, and J. G. Grossman, “Extraor¬ 
dinary sunlight absorption and one nanometer thick photo¬ 
volt aics using two-dimensional monolayer materials,” Nano 
letters, vol. 13, no. 8, pp. 3664-3670, 2013. 

O. Lopez-Sanchez, D. Lembke, M. Kayci, A. Radenovic, 
and A. Kis, “Ultrasensitive photodetectors based on mono- 
layer M 0 S 2 ,” Nature nanotechnology, vol. 8, no. 7, pp. 497- 
501, 2013. 

J. S. Ross, P. Klement, A. M. Jones, N. J. Ghimire, 
J. Yan, D. Mandrus, T. Taniguchi, K. Watanabe, K. Kita- 
mura, W. Yao, et al, “Electrically tunable excitonic light- 
emitting diodes based on monolayer WSe 2 pn junctions,” 
Nature nanotechnology, vol. 9, no. 4, pp. 268-272, 2014. 

A. Pospischil, M. M. Furchi, and T. Mueller, “Solar-energy 
conversion and light emission in an atomic monolayer pn 
diode,” Nature nanotechnology, vol. 9, no. 4, pp. 257-261, 

2014. 

G. Grosso and G. Parravicini, Solid State Physics. Elsevier 
Science, 2000. 

P. Gudazzo, G. Attaccalite, 1. V. Tokatly, and A. Ru¬ 
bio, “Strong charge-transfer excitonic effects and the bose- 
einstein exciton condensate in graphane,” Phys. Rev. Lett., 
vol. 104, p. 226804, Jun 2010. 

P. Gudazzo, 1. V. Tokatly, and A. Rubio, “Dielectric screen¬ 
ing in two-dimensional insulators: Implications for exci¬ 
tonic and impurity states in graphane,” Phys. Rev. B, 
vol. 84, p. 085406, Aug 2011. 

F. Hiiser, T. Olsen, and K. S. Thygesen, “How dielec¬ 
tric screening in two-dimensional crystals affects the con¬ 
vergence of excited-state calculations: Monolayer M 0 S 2 ,” 
Phys. Rev. B, vol. 88, p. 245309, Dec 2013. 


O. Pulci, P. Gori, M. Marsili, V. Garbuio, R. Del Sole, and 
F. Bechstedt, “Strong excitons in novel two-dimensional 
crystals: Silicane and germanane,” EPL (Europhysics Let¬ 
ters), vol. 98, no. 3, p. 37004, 2012. 

T. G. Berkelbach, M. S. Hybertsen, and D. R. Reich- 
man, “Theory of neutral and charged excitons in monolayer 
transition metal dichalcogenides,” Phys. Rev. B, vol. 88, 
p. 045318, Jul 2013. 

H. Terrones, F. Lopez-Urias, and M. Terrones, “Novel 
hetero-layered materials with tunable direct band gaps 
by sandwiching different metal disulfides and diselenides,” 
Scientific reports, vol. 3, 2013. 

L. Britnell, R. Ribeiro, A. Eckmann, R. Jalil, B. Belle, 
A. Mishchenko, Y.-J. Kim, R. Gorbachev, T. Georgiou, 
S. Morozov, et al., “Strong light-matter interactions in het¬ 
erostructures of atomically thin films,” Science, vol. 340, 
no. 6138, pp. 1311-1314, 2013. 

A. Geim and 1. Grigorieva, “Van der waals heterostruc¬ 
tures,” Nature, vol. 499, no. 7459, pp. 419-425, 2013. 

K. Andersen, S. Latini, and K. S. Thygesen, “Dielectric 
genome of van der waals heterostructures,” Nano letters, 
vol. 15, no. 7, pp. 4616-4621, 2015. 

S. L. Adler, “Quantum theory of the dielectric constant in 
real solids,” Physical Review, vol. 126, no. 2, p. 413, 1962. 

T. Gheiwchanchamnangij and W. R. L. Lambrecht, “Quasi¬ 
particle band structure calculation of monolayer, bilayer, 
and bulk M 0 S 2 ,” Phys. Rev. B, vol. 85, p. 205302, May 
2012. 

A. Molina-Sanchez, D. Sangalli, K. Hummer, A. Marini, 
and L. Wirtz, “Effect of spin-orbit interaction on the op¬ 
tical spectra of single-layer, double-layer, and bulk mos 2 ,” 
Phys. Rev. B, vol. 88, p. 045412, Jul 2013. 

“The dielectric building blocks and the qeh software can 
be downloaded from:.” https://cmr.fysik.dtu.dk/vdwh/ 
vdwh.html. Accessed: 2015-08-18. 

F. A. Rasmussen and K. S. Thygesen, “Gomputational 
2D materials database: Electronic structure of transition 
metal dichalcogenides and oxides,” The Journal of Physical 
Chemistry C, 2015. 

F. Hiiser, T. Olsen, and K. S. Thygesen, “Quasiparticle 
gw calculations for solids, molecules, and two-dimensional 
materials,” Phys. Rev. B, vol. 87, p. 235132, Jun 2013. 

G. Strinati, “Effects of dynamical screening on resonances 
at inner-shell thresholds in semiconductors,” Phys. Rev. B, 
vol. 29, pp. 5718-5726, May 1984. 

G. Onida, L. Reining, and A. Rubio, “Electronic ex¬ 
citations: density-functional versus many-body green’s- 
function approaches,” Rev. Mod. Phys., vol. 74, pp. 601- 
659, Jun 2002. 

G. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and 
A. Rubio, “Exact coulomb cutoff technique for supercell 
calculations,” Phys. Rev. B, vol. 73, p. 205119, May 2006. 
J. Enkovaara, G. Rostgaard, J. J. Mortensen, J. Ghen, 

M. Dulak, L. Ferrighi, J. Gavnholt, G. Glinsvad, 
V. Haikola, H. Hansen, et al, “Electronic structure cal¬ 
culations with gpaw: a real-space implementation of the 
projector augmented-wave method,” Journal of Physics: 
Condensed Matter, vol. 22, no. 25, p. 253202, 2010. 

“The gpaw code is available as part of the campos soft¬ 
ware:.” https://wiki.fysik.dtu.dk/gpaw/. Accessed: 
2010-09-30. 



14 


J. Yan, K. W. Jacobsen, and K. S. Thygesen, “Opti¬ 
cal properties of bulk semiconductors and graphene/boron 
nitride: The bethe-salpeter equation with derivative 

discontinuity-corrected density functional energies,” Phys. 
Rev. B, voL 86, p. 045208, Jul 2012. 

A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, 
Y. Li, O. B. Aslan, D. R. Reichman, M. S. Hybertsen, 
and T. F. Heinz, “Exciton binding energy and nonhydro- 
genic rydberg series in monolayer WS 2 ,” Phys. Rev. Lett.^ 


vol. 113, p. 076802, Aug 2014. 

H. Peelaers and C. G. Van de Walle, “Effects of strain on 
band structure and effective masses in M 0 S 2 ,” Phys. Rev. 
B, vol. 86, p. 241401, Dec 2012. 

H.-P. Komsa and A. V. Krasheninnikov, “Effects of con¬ 
finement and environment on the electronic structure and 
exciton binding energy of M 0 S 2 from first principles,” Phys. 
Rev. B, vol. 86, p. 241201, Dec 2012. 



